Raspberry Pi 5 + Mathematica: 
A STEM Creativity Powerhouse

Em. prof. Bart ter Haar Romeny
Eindhoven, the Netherlands

#WolframTechConf

ABSTRACT

Discover how the Raspberry Pi 5 turns into a full STEM lab with free Mathematica. From real-time sensor dashboards to image processing and stunning math art, this guide empowers students, teachers, and parents to explore, invent, and visualize science and engineering like never before.

The Raspberry Pi 5 - a STEM student’s dream

STEM: Science, Technology, Engineering, Mathematics.

Raspberry Pi is a series of small single-board computers (SBCs) developed in the UK by the Raspberry Pi Foundation in collaboration with Broadcom.

Originally created to help teach computer science in schools, but it is now used in areas such as industrial automation, robotics, home automation, IoT devices, and hobbyist projects.

Mathematica is free on every Raspberry Pi

Mathematica is IDEAL for STEM students and teachers, but the old models (RPI3, RPI4) were too weak to run Mathematica

Not a single Raspberry Pi project website / vendor / community mentions the use of Mathematica

But now: the Raspberry Pi 5 is powerful enough, and Mathematica 14.2 runs amazingly well

Every function is included (incl. all math, parallel processing, compilation, camera, full interactivity, 3D visualizations etc.), except CUDA support and some video functions.

The RPI-5 features:
- 2.4 GHz quad-core Cortex-A76 CPU
- a VideoCore VII GPU
- PCIe support, and 2/4/8/16 GB RAM

8 or 16 GB (best) of RAM ($ 100 / $ 160 on Amazon) is a required minimum to run Mathematica

Save money: with the free NoMachine software a bare bone RPI5 can be remotely controlled by any (even old) laptop, desktop or tablet.

How to prepare your RPI5: https://www.21stcenturycoding.net/prepare/

Installation is with two clicks!

Ubiquitous and popular: 
1.9 million RPI-5 units sold in 2024

Interact with an amazing array of (cheap) sensors, motors, displays, joysticks, cameras, LEDs, speakers, switches, pumps etc. through GPIO pins, Bluetooth, USB ports and WiFi

($ 15)

Most users struggle with cumbersome software (often outdated and old-fashioned, with endless searches on many user groups)

Integration of 
- low-level Python/C-C++ for interface control and 
- high-level Mathematica for analysis and visualization
gives a super-boost to projects

shape=Import["https://www.neuromath.net/images/RPI5.obj"];
Graphics3D[shape]

Mathematica is cheap and mature

Mathematica is free on every Raspberry Pi (URL).

Mathematica is free for most university students (URL)

Mathematica is cheaper than Netflix for all STEM students ($ 7,35/month, URL) on any other computer

Current LLMs today typically learn Python, C++ and Java. It would be timely to add the WL code of millions of new Mathematica users!

Let me amaze you with this presentation

All applications of this presentation are executed
on a Raspberry Pi 5 with 16 GB and 256 GB NVMe.

Detailed instructions how to prepare your Raspberry Pi 5

Installation RPI5 Operating System, installation Mathematica, installation remote viewing software, installation MVMe disk, running benchmarks:
Visit www.21stcenturycoding.net/prepare.
For a 1 hour ‘Gentle Introduction to Mathematica’ and more applications see: www.youtube.com/@21stcenturycoding and www.21stcenturycoding.net .

About me

I always was a STEM student. I used Mathematica from version 1, almost daily. Now 73, retired professor in medical image analysis at Eindhoven University of Technology. I did all my teaching and research in Mathematica. In 2013 I received the Mathematica Innovator Award.

I told my students: If your work is your hobby, you never have to work again. 
Mathematica was and is my hobby, it made my career.

$Version

SetDirectory[NotebookDirectory[]];

Introduction

I selected here a few applications that show the versatility of the Wolfram Language specifically for the RPI5, with explanation, in computational essay style, for STEM students.

No ‘wow, I can write this with an LLM’, but focused on understanding and learning in depth. Wolfram’s material stack is huge: https://www.wolfram.com/language/?source=nav#learning-resources

Sensors

Temperature, pressure and humidity with gauges, LCD display and a 1$ sensor

The sensor BME280 measure temperature, pressure and humidity. 
The LCD1602 display has 4 rows of 20 characters.
We read the values from the sensor and write the characters to the display with the same standard I2C connection of the RPI5.

Requirements: 
Raspberry Pi 5 with 8 or (preferably) 16 GB RAM
BME280 sensor (e.g. Amazon)
4 row - 20 chars LCD display (e.g. Amazon) 
A breadboard and connection cables.

The I2C interface

The I2C communication protocol is widely used for simple communication between processors and chips. It only requires two lines (SDA for data and SCL for the clock), next to 2 lines for  power. Learn more about I2C on Wikipedia: "https://en.wikipedia.org/wiki/I%C2%B2C"

To use the I2C, it needs to be enabled. Select on the top bar menu:
RPI logo -> Preferences -> Raspberry Pi Configuration -> Interfaces -> turn I2C on.

Connect the BME280 sensor, by 4 wires, to the following physical pin numbers:
Vcc -> PIN 1 (NB: the BME280 sensor is a 3.3 Volt device!)
Gnd -> PIN 9
SDA -> PIN 3
SCL -> PIN 5

The next cell is an ExternalLanguage cell. Enter is by typing > as the first character, and select with your mouse in the cell bullet on the left side of the cell which type of language you want. Here we select the (Linux) Shell.

This shows the pinout of your Raspberry Pi 5:

pinout

Then, connect the LCD display, by 4 wires:
Vcc -> PIN 2 (The LCD1602 is a 5V device.)
Gnd -> PIN 9
SDA -> PIN 3
SCL -> PIN 5

After adding the LCD1602 display to the same I2C (yes, you can add multiple sensors to a single I2C channel), we see that it has address 27 (this was specified in the manual ...) and the I2C address of the BME280 sensor is 76:

We can check :

i2cdetect -y 1

There are numerous webpages describing how to use sensors and the Raspberry Pi (e.g. https://learn.adafruit.com/).

Here I show how particularly easy things are in Mathematica, exploiting pieces of Python code to interact with the external devices.

The Python session for the hardware

We use Python to read the sensor and write the display, and use Mathematica for the main program, and possible fancy interfaces and data analysis.

session=StartExternalSession["Python"]

Note that a ‘virtual environment’ Venv is created automatically by Mathematica. You don’t have to worry about directories or paths, it is all taken care of by Mathematica.

Install required Python packages

Installation of the required Python packages is a one-liner in Mathematica. Don’t forget to bundle multiple packages in a List with curly brackets, as is the standard list format in Mathematica.

This loads the required drivers for the BME280 sensor and the LCD display:

ResourceFunction["PythonPackageInstall"][{"RPI.BME280","RPLCD"}]

Read the sensor

This Python code, run in the ExternalEvaluate[ ] command in the just defined session initializes the BME280 sensor (packages are imported, the address is set, the bus is initialized):

ExternalEvaluate[session,"
import time
import smbus2
import bme280

# BME280 sensor address (default address)
address = 0x76

# Initialize I2C bus
bus = smbus2.SMBus(1)

# Load calibration parameters
calibration_params = bme280.load_calibration_params(bus, address)
"]

We then define a Mathematica function in Python, which we can call in Mathematica.
The Python code ‘def read()’ defines a function read() in Python. The function has no input argument(s), and returns a triplet of measurement values :

readSensor:=ExternalEvaluate[session,"
# Read sensor data
data = bme280.sample(bus, address, calibration_params)
def read():
    temperature_celsius = data.temperature
    pressure = data.pressure
    humidity = data.humidity
    return (temperature_celsius, pressure, humidity)
"]

We defined the function readSensor in Mathematica with := (‘delayed assignment).

It is an ExternalFunction. You can check by asking information about readSensor.

Information[readSensor]

Now call the function readSensor[ ] in Mathematica (the square brackets are obligatory, as it is a function call):

{temp,press,hum}=readSensor[]

Quite nice, we get the 3 readings with a Mathematica call. 
This method can now be applied to any I2C sensor.

Write the display

This initializes the LCD display:

ExternalEvaluate[session,"
from RPLCD.i2c import CharLCD
from time import sleep

# Initialize 20x4 LCD on I2C bus 1 at address 0x27
lcd = CharLCD(i2c_expander='PCF8574', address=0x27, port=1,
              cols=20, rows=4, charmap='A02',
              auto_linebreaks=True)
"]

And here is the function to write the date, temperature, pressure and humidity on each line.
This function takes 4 arguments:

writeLCD1602:=ExternalEvaluate[session,"
def write(dateString,temp,press,hum):
#    lcd.clear()
    # Set cursor position for date
    lcd.cursor_pos = (0, 0)
    lcd.write_string(dateString)

    # Set cursor position for temperature
    lcd.cursor_pos = (1, 1)
    temp_str = f'Temp :   {temp:.2f} C'
    lcd.write_string(temp_str)

    # Set cursor position for pressure
    lcd.cursor_pos = (2, 1)
    press_str = f'Press: {press:.2f} hPa'
    lcd.write_string(press_str)

    # Set cursor position for humidity
    lcd.cursor_pos = (3, 1)
    hum_str = f'Humid:   {hum:.2f} %'
    lcd.write_string(hum_str)
"]

Create the date string:

dateString=DateString[{"DayNameShort"," ","Year","-","Month","-","Day"," ","Hour",":","Second"}]

and test the function:

writeLCD1602[dateString,temp,press,hum]

Dynamically update every second

Finally, we dynamically read the sensor and update the LCD display and a GUI with gauges every second:

DynamicModule[{},
Dynamic[dateString=DateString[{"DayNameShort"," ","Year","-","Month","-","Day"," ","Hour",":","Minute"}];
{temp,press,hum}=readSensor[];
writeLCD1602[dateString,temp,press,hum];
Framed[Row[{ClockGauge[],ThermometerGauge[temp,{15,35}],Column[{AngularGauge[press,{1000,1050},ScaleOrigin->{π,0}],AngularGauge[hum,{0,100},ScaleOrigin->{π,0},PlotTheme->"Scientific"]}],Column[{
dateString,
" Temp :   "<>ToString[NumberForm[temp,{∞,2}]]<>" C",
" Press: "<>ToString[NumberForm[press,{∞,2}]]<>" hPa",
" Humid:   "<>ToString[NumberForm[hum,{∞,2}]]<>" %" 
}]}]],UpdateInterval->1]]

Take the BME280 sensor between your fingers, and see the temperature change interactively.

In Mathematica (Wolfram Language) it is really easy to make fancy user interfaces in a few lines.

Note that the DynamicModule (as does any Manipulate function) pauses when the output window is scrolled out of sight. This then also stops the updating of the LCD1602 display.
Stop the system by deleting the above cell with the gauges.

Blush, baby, blush

We can now interactively let a girl blush.

The sensor readings are used to make the colour of -only those pixels in the face that represent skin- more reddish. This more complex image analysis is easy in Mathematica, and real-time.

Face selection and skin redness

To work with an anonymous face, the website in the next statement comes handy. Each time it is called, a new random face is generated, 1024 x 1024 pixels.

face=Import["https://www.thispersondoesnotexist.com","Image"]

For the demo notebook, we use the face below, but you can of course take any face.

face=;

Check the image dimensions:

ImageDimensions[face]

A neural network in Mathematica can recognize some facial features:

FacialFeatures[face]

We now go find the RGB color of the skin. First we put the image on screen:

face

Put your mouse pointer on the location where you want to measure the skin color, e.g. on the center of the cheek. Press the right mouse button, and select ‘Get Coordinates’. Click the left mouse key, and Ctrl-c to copy these coordinates.

Paste these coordinates into the highlighted section of the next statement:

ColorConvert[PixelValue[face,{702.6709349092998,366.5915922943212}],OptionValue[Options[face],ColorSpace]->ImageColorSpace[face]]

It gives you the color, as a colored pixel. If you really need the RGB values, use:

//FullForm

You can now copy-paste the coloured pixel into the next statement, which finds all pixels with a color near the specified color, with a ‘color-distance’ in the 3D RGB space of 0.01, and recolours them to Red.

reddish=ImageRecolor[face,ColorsNear[,0.01]->Red]

Nice. Only the skin pixels are modified, quite a sophisticated result.

Blush

The girl blushes by adding more reddish when the temperature of the sensor rises.
Therefore we first need to calibrate the system, as it is dependent on your local room temperature (the starting temperature), and the final temperature (e.g. by putting your fingers on both side of the sensor, or to blow warm air over it.

With Dynamic we get a continuous reading of the temperature. Write down the start and end temperature, e.g. 19.1 and 19.8 degrees.

The overlay opacity should change from 0.1 to 0.5, corresponding with the begin and end temperature.

Then fit a straight line (use in Fit only the constant term 1, and the linear term x) between the points {18,0.1} and {22,0.5}, which gives us the offset and slope:

Fit[{{18,0.1},{22,0.5}},{1,x},x]

Enter this in the following statement:

Dynamic[Show[ImageCompose[face,{reddish,-1.70+0.1 First@readSensor[]}],ImageSize->300]]
Dynamic[First@readSensor[]]

She will blush when you slightly heat the BME280 sensor with your fingers.

YouTube: https://youtu.be/jXFhIfibOHY

Compile

MathArt - swimming jellyfish in an aquarium

A tribute to Tom Wickham - Jones, (Director of Kernel Technology at Wolfram Research)

Source: https://community.wolfram.com/groups/-/m/t/3320126, 14-11-2024

This example is taken from the Wolfram Community posts, where Vitaliy Kaurov, Director of Engagement at Wolfram Research, discovered a nice set of mathematical art implementations by @yuruyurau written in JavaScript. He took the example to show how ChatGPT could transfer the Java code into Wolfram Language.

As the Raspberry Pi 5 is a tiny computer, we go make some adjustments to make it faster and suitable for our RPI5.

But above that, this notebook will learn you a nice set of programming paradigms.

Here is the original mathematical function:

a[x_, y_, t_] := Module[{k, e, d, q, c},
  k = x/8 - 25;
  e = y/8 - 25;
  d = (k^2 + e^2)/99;
  q = x/3 + k/Cos[5y] Sin[d^2 - t]/2;
  c = d/2 - t/8;
  {
    q*Sin[c] + e*Sin[d + k - t] + 200,
    (q + y/8 + d*9)*Cos[c] + 200
  }
]

In the code the math formula is in a nested form, and it is difficult to see what happens. The resulting formula is quite complicated. So let’s view it in the traditional (text book) form:

a[x,y,t]//TraditionalForm

Pretty amazing indeed. It is a function that gives as output a pair of values at coordinates {x,y} at time t. The input consists of three variables, x and y for the position and t for the time. We go make a canvas of 400 x 400 positions (pixels), and we go scan this canvas for all positions 100 ≤ x ≤ 299 in steps of 1, and 99 ≤ y ≤ 299 in steps of 5.

Let us inspect what the function a[x,y,t] actually gives. Here is the value of the first time point at t=0 at the starting position {x=100, y=99}:

a[100,99,0]

Oeps. We immediately see the problem: as Mathematica works with symbolic expressions when exact numbers (like integers) are involved, this is an expensive operation, both in memory use and time. The situation is sped up when we take real numbers. This is simply done by converting an integer by adding a numerical point after it, so t=0 becomes t=0. or by using the numerical value function N[0].

a[100,99,0.]

We only have to specify one real to force Mathematica to process the formula numerically:

a[100.,99,0]

Here is a timing experiment, to calculate the complete jellyfish image at 10 time moments:

time=First@Timing[pts=Table[a[x,y,t],{x,99,299,5},{y,100,299},{t,0,9}];]

And with reals:

timeNumerical=First@Timing[pts=Table[a[x,y,t],{x,99,299,5},{y,100,299},{t,0.,9.}];]

A bit of gain.

Can we do more by compiling the function? Let us do the exercise, and learn how compilation works.

When a function is compiled, the code is made into the lowest possible machine code, and with pre-specified number types. See for more details the Documentation.

The Compile function is easy to build. For example, f = Compile[{{x,_Real},{y,_Real}},x^2+4 y^2] gives the function f[x,y] that calculates (x^2)+4(y^2). We need to give the type of the internal representation of the variables (both x and y are Real), and do that first within the Compile function brackets.

We convert the formula to InputForm so it can be used as code in Compile:

a[x,y,t]//InputForm

We copy it into the Compile body, also changing the non-exponent integer numbers to Reals (by adding the numerical dot):

aCompiled=Compile[{{x,_Real},{y,_Real},{t,_Real}},
{200.-(-25.+y/8.)*Sin[25.+t-x/8.+(-(-25.+x/8.)^2-(-25.+y/8.)^2)/99.]-Sin[t/8.+(-(-25.+x/8.)^2-(-25.+y/8.)^2)/198.]*(x/3.-((-25.+x/8.)*Sec[5.*y]*Sin[t-((-25.+x/8.)^2+(-25.+y/8.)^2)^2/9801.])/2.),200.+Cos[t/8.+(-(-25.+x/8.)^2-(-25.+y/8.)^2)/198.]*(x/3.+((-25.+x/8.)^2+(-25.+y/8.)^2)/11.+y/8.-((-25.+x/8.)*Sec[5.*y]*Sin[t-((-25.+x/8.)^2+(-25.+y/8.)^2)^2/9801.])/2)}]

And the explanation: it is a function of three arguments, and we need to specify the type of both the arguments and the internally used variables. The typing for the arguments is done at the very beginning: a list (specified by the curly brackets) of lists: {{x,_Real},{y,_Real},{t,_Real}}. The underscore _ means ‘any value’, and _Real means ‘any value, typed as Real’. Then the main body specifies the formula. The typing of internally used variables is done after the main code (see the documentation of Compile), but we don’t have them here. They would have been specified after a comma after the main code, in the same format as the argument types.

The proof of the pudding is in the eating, so here is the timing test, for 10 frames:

timeCompiled=First@Timing[pts=Table[aCompiled[x,y,t],{x,99,299,5},{y,100,299},{t,0,9}];]

How much faster?

(timeNumerical/timeCompiled)

Quite nice!

A remark: why not put the original code in the Compile function? Well, as this code is quite nested, each time a complicated calculation has to be done for each point. The strength of Mathematica is that it can do this expansion of the nesting once, symbolically.

Here is the timing result if we put the nested code into the Compile function:

aCompiledNested=Compile[{{x,_Real},{y,_Real},{t,_Real}},
k = x/8. - 25.;
  e = y/8. - 25.;
  d = (k^2. + e^2.)/99.;
  q = x/3. + k/Cos[5.y]Sin[d^2. - t]/2.;
  c = d/2. - t/8.;
  out={
    q*Sin[c] + e*Sin[d + k - t] + 200.,
    (q + y/8. + d*9.)*Cos[c] + 200.
  };
out,
{{k,_Real},{e,_Real},{d,_Real},{q,_Real},{c,_Real},{out,_Real,1}}]

timeCompiledNested=First@Timing[pts=Table[aCompiledNested[x,y,t],{x,99,299,5},{y,100,299},{t,0,10}];]

(timeNumerical/timeCompiledNested)

So that could hardly beat Mathematica’s own implementation ...

The drawing function:

Take a background image, e.g. an aquarium, preferably with dark appearance, so the white dots of the jellyfish stand out nicely:

bg=ImageResize[ImageTake[Import["https://thumbs.dreamstime.com/b/saltwater-full-corals-mixed-aquarium-one-most-beautiful-addition-home-coral-reef-aquarium-tank-scene-fishes-118136040.jpg"],{All},{254,800}],{400,400}]

drawc[t_]:=Module[{pointsList},
pointsList=Flatten[Table[aCompiled[x,y,t],{y,99,299,5},{x,100,299}],1];
Graphics[{{Opacity[0.5],White,PointSize[.002],Point[pointsList]}},
Prolog -> Inset[Graphics[bg],{0,0},{0,0},420],
AspectRatio->1,PlotRange->{{0,400},{0,400}}]]

Here is the interactive motion of the jellyfish:

Manipulate[drawc[t/10.],{t,0,160π,1}]

Click the + to the right of the slider to show the controls, and hit the play button.

Geogebra with the Wolfram Language - write your own apps easily

Circle from 3 points by perpendicular bisectors

Geogebra is a computer algebra system, written in JavaScript and HML5.

According to German Wikipedia, GeoGebra had about 30 million registered users across 190 countries by 2017, increasing to around 100 million users in 195 countries by December 2021.

On Android, GeoGebra apps have amassed approximately 30 million total installs across its suite.

To find the equation for a circle, going through three points, we use the fact that two bisectors (the line halfway and normal to the line connecting a pair of points) intersect in the center point of the circle.

Carefully try to understand the code below, and play with the resulting animated figure.

You can turn the bisectors on and off in the display.

Manipulate[
(*Define the 3 points*)
{p,q,r}=locators;
(*Calculate the center points of PQ and QR*)
midPQ=(p+q)/2;
midQR=(q+r)/2;
(*Calculate the slopes of the lines*)
slopePQ=If[q[[1]]==p[[1]],∞,((q[[2]]-p[[2]])/(q[[1]]-p[[1]]))];slopeQR=If[r[[1]]==q[[1]],∞,((r[[2]]-q[[2]])/(r[[1]]-q[[1]]))];
(*Define the perpendicular bisectors*)
perpSlopePQ=If[slopePQ===∞,0,If[slopePQ==0,∞,-(1/slopePQ)]];perpSlopeQR=If[slopeQR===∞,0,If[slopeQR==0,∞,-(1/slopeQR)]];bisectorPQ=If[perpSlopePQ===∞,x==midPQ[[1]],y-midPQ[[2]]==perpSlopePQ (x-midPQ[[1]])];bisectorQR=If[perpSlopeQR===∞,x==midQR[[1]],y-midQR[[2]]==perpSlopeQR (x-midQR[[1]])];center=Solve[{bisectorPQ,bisectorQR},{x,y}];
If[center==={},Return[Graphics[Text["Points are collinear"]]]];
centerx=x/.center[[1]];
centery=y/.center[[1]];
(*Calculate the radius*)
radius=(Sqrt[(((centerx-p[[1]]))^2)+(((centery-p[[2]]))^2)]);
(*Draw the circle and points*)
Show[{Graphics[{
Line[{p,q}],
Line[{q,r}],
If[wl,{Texture[],
{Disk[{centerx,centery},radius],Circle[{centerx,centery},radius]}},
Circle[{centerx,centery},radius]],
Line[{p,q}],
Line[{q,r}],
Red,PointSize[Large],Point[p],Point[q],Point[r],
Green,Point[{x,y}/.center]
},PlotRange->{{-10,10},{-10,10}},AspectRatio->1],If[showbisector1,Plot[y/.Solve[bisectorPQ,y],{x,-10,10}],Graphics[]],If[showbisector2,Plot[y/.Solve[bisectorQR,y],{x,-10,10}],Graphics[]]}],
(*Manipulate locators*)
{{locators,{{-3,-4},{3,-3},{4,5}}},Locator},
{{showbisector1,False,"Show perpendicular bisector 1"},{True,False}},{{showbisector2,False,"Show perpendicular bisector 2"},{True,False}},
{{wl,False,"Wolfram Logo"},{True,False}},
FrameLabel->{"x","",Style["Circle through 3 points",35]}
]

How to plot a perpendicular bisector?

bisectorPQ

We want y as function of x, so we solve for y:

Solve[bisectorPQ,y]

Plot[y/.Solve[bisectorPQ,y],{x,-10,10},PlotRange->{{-10,10},{-10,10}},AspectRatio->1]

The world’s easiest coding with Neural Nets -

Interactive object detection with a camera and Yolo5

The 165 neural networks in the Wolfram Neural Network Repository are an amazing toolkit for neural net experimentation.

NetModel[]

How many networks are currently available?

Length[NetModel[]]

Let’s do object recognition (detection -> localization, and classification -> naming) with a regular USB webcam hooked up to our Raspberry Pi 5. Connect it to one of the two blue USB-3 connectors, which are faster.

It will not be real-time, but, as you will see, we get about 1 frame/second on the CPU of the Raspberry Pi 5. Quite nice.

No software and hardware hassle, no expensive AI hat, no complicated camera drivers: it just works!
With just a few lines of code.

net=NetModel["YOLO V5 Trained on MS-COCO Data"]

We exploit the code given in the worked-out example on the repository website for this network:

nonMaximumSuppression=ResourceFunction["NonMaximumSuppression"];

labels={"person","bicycle","car","motorcycle","airplane","bus","train","truck","boat","traffic light","fire hydrant","stop sign","parking meter","bench","bird","cat","dog","horse","sheep","cow","elephant","bear","zebra","giraffe","backpack","umbrella","handbag","tie","suitcase","frisbee","skis","snowboard","sports ball","kite","baseball bat","baseball glove","skateboard","surfboard","tennis racket","bottle","wine glass","cup","fork","knife","spoon","bowl","banana","apple","sandwich","orange","broccoli","carrot","hot dog","pizza","donut","cake","chair","couch","potted plant","bed","dining table","toilet","tv","laptop","mouse","remote","keyboard","cell phone","microwave","oven","toaster","sink","refrigerator","book","clock","vase","scissors","teddy bear","hair drier","toothbrush"};

netevaluate[model_,img_,detectionThreshold_:.5,overlapThreshold_:.5]:=Module[{imgSize,classes,coords,obj,scores,bestClass,probable,probableClasses,probableScores,probableBoxes,h,w,max,scale,padding,nms,finals},imgSize=Last@NetExtract[model,{"Input","Output"}];
{classes,coords,obj}=Values@model[img];
(*each class probability is rescaled with the box objectness*)scores=classes*obj;
bestClass=Last@*Ordering/@scores;
(*filter by probability*)(*very small probability are thresholded*)probable=UnitStep[obj-detectionThreshold];{probableClasses,probableBoxes,probableScores}=Map[Pick[#,probable,1]&,{labels[[bestClass]],coords,obj}];
If[Length[probableBoxes]==0,Return[{}]];
(*transform coordinates into rectangular boxes*){w,h}=ImageDimensions[img];
max=Max[{w,h}];
scale=max/imgSize;
padding=imgSize*(1-{w,h}/max)/2;
probableBoxes=Apply[Rectangle[scale*({#1-#3/2,imgSize-#2-#4/2}-padding),scale*({#1+#3/2,imgSize-#2+#4/2}-padding)]&,probableBoxes,1];
(*gather the boxes of the same class and perform non-max suppression*)nms=nonMaximumSuppression[probableBoxes->probableScores,"Index"];
finals=Transpose[{probableBoxes,probableClasses,probableScores}];
Part[finals,nms]];

Mathematica’s function CurrentImage[] reads an image from the USB webcam (or any other connected camera)

testImage=CurrentImage[];
detection=netevaluate[net,testImage]
HighlightImage[testImage,GroupBy[detection[[All,1;;2]],Last->First]]

DynamicModule[{currentImage,detection},
  Manipulate[
    currentImage=CurrentImage[];
    detection=netevaluate[net,currentImage];
    HighlightImage[currentImage,GroupBy[detection[[All,1 ;; 2]],Last->First],
      ImageSize->Medium
    ],
    {{Automatic,"Update"},None,"Update"},
    ControlType->Trigger,
    TrackedSymbols:>{currentImage}
  ]
]

Quite impressive, one full frame detection and classification per second, on your tiny Raspberry Pi 5!

YouTube video: https://youtu.be/zAaxj2KOy1o

Complex Neural Net examples

Mathematica has a professional neural network environment. Everything is possible and easy. Here a complex net is constructed by telling how the connect the outputs of the numbered layers to the inputs of other layers, it is intuitive:

NetGraph[{BatchNormalizationLayer[],Tanh,LogisticSigmoid,Tanh,TotalLayer[],TotalLayer[],TotalLayer[],CatenateLayer[],LinearLayer[50]},{1->2,1->3,1->4,2->5,3->5,3->6,4->6,2->7,4->7,5->8,6->8,7->8,8->9},"Input"->20]

Expand the graph by clicking on the +. Clicking on any network element gives you the full information of what goes in and out, and what is the size of the data there.

And how about this one?

See how fast your RPI5 handles this exquisite task (again click on the + to visualize the network and use the scrollbar on the bottom of the screen)):

NetInitialize[NetGraph[{BatchNormalizationLayer[],Tanh,LogisticSigmoid,Tanh,TotalLayer[],TotalLayer[],TotalLayer[],CatenateLayer[],LinearLayer[50],BatchNormalizationLayer[],Tanh,LogisticSigmoid,Tanh,TotalLayer[],TotalLayer[],TotalLayer[],CatenateLayer[],LinearLayer[50],DropoutLayer[],DropoutLayer[],TotalLayer[],LogisticSigmoid,BatchNormalizationLayer[],Tanh,LogisticSigmoid,Tanh,TotalLayer[],TotalLayer[],TotalLayer[],CatenateLayer[],LinearLayer[50],DropoutLayer[],TotalLayer[],TotalLayer[],TotalLayer[],TotalLayer[],TotalLayer[],TotalLayer[],LinearLayer[50],LinearLayer[4]},{1->2,1->3,1->4,2->5,3->5,3->6,4->6,2->7,4->7,5->8,6->8,7->8,8->9,10->11,10->12,10->13,11->14,12->14,11->15,13->15,13->16,12->16,16->17,15->17,14->17,17->18,18->19,9->20,20->21,19->21,21->22,23->24,23->25,23->26,24->27,25->27,24->28,25->29,26->28,26->29,27->30,28->30,29->30,30->31,31->32,32->21,30->33,8->33,8->34,17->34,30->35,17->35,33->36,34->36,34->37,35->37,37->38,36->38,38->39,39->21,22->40},"Input"->85]]

Visualize all internal layers of a neural net interactively

If we want to know what is internally represented in the hidden layers of the network, we can easily visualize these layers in Mathematica.
Let’s take as example this neural net, and open the resulting NetChain by pressing the + :

NetModel["Wolfram ImageIdentify Net V1"]

(CopyTag[NetChain[]])

You see that there are many layers, and each has a name: the first layer is called conv_1, the second layer is called bn_1 etc. The input image is an arbitrary image of a cat (see how nice it is to have images in your code ☺).

The following command shows the list of internal images, per layer, in a Manipulate function. Click on the + sign to the right of the slider to open the settings, so you can easily do the inspection layer by layer:

Manipulate[GraphicsGrid[Partition[Image/@Take[NetModel["Wolfram ImageIdentify Net V1"],1][],n]],{{n,4,"Layer number"},1,30,1}]

The Machine Learning book by Etienne Bernard

This is the best tutorial book on machine learning and neural nets I have seen. Fully interactive, covers a wide range, written by one of the main authors of Neural Nets in the Wolfram Language. All notebook chapters can be downloaded for free.

Import["https://www.wolfram-media.com/products/img/introduction-to-machine-learning.png"]

Computer vision

Smile, Mona Lisa, smile!

Here I show that you can have fast interactivity on images, which means thousands of pixels, on the RPI-5. We go interactively deform the face of Mona Lisa, by moving points at specific positions, and interpolate smoothly the new deformed area around the shifted point. The method for such a radial interpolation is called ‘radial basis interpolation’, using the ‘radial basis transform’. The code for this transform is given at the end, as initialization cell. Gratitude to Markus van Almsick for the code, who worked for 8 years in my Biomedical Image Analysis group at Eindhoven University of Technology, as PhD student and as postdoc, and later worked for many years as image processing expert at Wolfram research.

This is a DyamicModule, a routine that continuously runs. The ‘sources’ are the locations of the points on her face. These are used in the Manipulate function to move interactively with the Manipulate function ‘Locator’. They are placed as Points, in Red.

This notebook shows that even on a small device as the Raspberry Pi 5 interactive deformation of images is possible.

DynamicModule[
{img=,
sources ={{190,100},{155,95},{131,100},{147,126},{206,194},{125,194}}},
Manipulate[
Magnify[
ImageTransformation[
img,
RadialBasisTransform[targets->sources,"ThinPlateSpline"],
Padding-> Gray,DataRange-> Full
],
2
],
{{targets,sources},Locator,Appearance->{Graphics[{Red, Point[{0,0}]}]}}
]
]

Cool visualizations

Make the famous 3Blue1Brown visualizations with a few lines of code

The interactive graphics in the YouTube tutorials of 3Blue1Brown (Gary Grant, www.youtube.com/@3blue1brown) are quite impressive. 
His explanations are outstanding, and his lectures have millions of views.

Here is an example of his graphics:

‘Constrained optimization introduction’: https://www.youtube.com/watch?v=vwUV2IDLP8Q

He uses a proprietary package ‘manim’ for his visualizations. These are not easily reproduced, needs an endless and complicated installation, and the code is difficult to read for a non-expert.

Here I show how easy it is to make such complex and interactive graphics in Mathematica.
Mathematica can even much more ...

Note, by the way, that the writing of mathematical formulas is done by hand, something virtually all YouTube authors still do ...

After reading this notebook you will appreciate the notion of ‘computational essay’.

Basic writing style in Mathematica

A short rehearsal:
Brackets: ( ) for grouping, { } for lists, [ ] for functions.
Mathematica symbols: Ctrl-^ for superscript, EscpEsc for π, etc.

To save memory, we do not keep the numbered outputs:

$HistoryLength=0;

The 3D graphics

The first function is f(x,y)=(x^2)y, of which we we view the ContourPlot3D:
(the function f-axis is the z-axis)

ContourPlot3D[(x^2) y==z,{x,-2,2},{y,-2,2},{z,-2,2}]

Rotate every 3D plot with your mouse, and keep the Ctrl key pressed to zoom, and keep the ⇧ key pressed to pan. 
Check the context-sensitive menu with a right mouse click.

We don’t want the Mesh, so we add Options to the ContourPlot3D function:

ContourPlot3D[{(x^2) y==z},{x,-2,2},{y,-2,2},{z,-2,2},Mesh->False]

Use the color picker from the ‘Insert’ menu:

ContourPlot3D[{(x^2) y==z},{x,-2,2},{y,-2,2},{z,-2,2},Mesh->False,ContourStyle->RGBColor[0.501961, 1, 1]]

The contour surface also needs Specularity:

ContourPlot3D[{(x^2) y==z},{x,-2,2},{y,-2,2},{z,-2,2},Mesh->False,ContourStyle->{RGBColor[0.501961, 1, 1],Specularity[1,100]}]

Check how Lighting is done in the Documentation.

ContourPlot3D[{(x^2) y==z},{x,-2,2},{y,-2,2},{z,-2,2},Mesh->False,ContourStyle->{RGBColor[0.501961, 1, 1],Specularity[1,100]},Lighting->"Accent"]

Rotate the plot to appreciate the specularities.

We add the cylinder (x^2)+(y^2)==1:

ContourPlot3D[{(x^2) y==z,(x^2)+(y^2)==1},{x,-2,2},{y,-2,2},{z,-2,2},Mesh->False,ContourStyle->{RGBColor[0.501961, 1, 1],Specularity[1,100]},Lighting->"Accent"]

Make the cylinder Yellow and Opacity[0.4] :

ContourPlot3D[{(x^2) y==c,(x^2)+(y^2)==1},{x,-2,2},{y,-2,2},{c,-2,2.},Mesh->False,ContourStyle->{{RGBColor[0.501961, 1, 1],Specularity[1,100]},{Yellow,Opacity[0.4]}},Lighting->"Accent"]

The BoundaryStyle is easy to understand:

ContourPlot3D[{(x^2) y==z,(x^2)+(y^2)==1},{x,-2,2},{y,-2,2},{z,-2,2},Mesh->False,ContourStyle->{{RGBColor[0.501961, 1, 1],Specularity[1,100]},{Yellow,Opacity[0.4]}},Lighting->"Accent",BoundaryStyle->{1->None,2->None,{1,2}->Red}]

The line is somewhat coarse, so we need more PlotPoints. 
And the boundary is not thick enough. We just add the style Thick.

ContourPlot3D[{(x^2) y==z,(x^2)+(y^2)==1},{x,-2,2},{y,-2,2},{z,-2,2},Mesh->False,PlotPoints->20,ContourStyle->{{RGBColor[0.501961, 1, 1],Specularity[1,100]},{Yellow,Opacity[0.4]}},Lighting->"Accent",BoundaryStyle->{1->None,2->None,{1,2}->{Red,Thick}}]

The 3 coordinate arrows are made as follows:

arrows=Graphics3D[{Arrowheads[.05],Arrow[Tube[{{0,0,0},{0,0,1.5}}]],
Arrow[Tube[{{0,0,0},{0,1.5,0}}]],
Arrow[Tube[{{0,0,0},{1.5,0,0}}]]}]

Any any set of graphics can be combined into the same coordinate system with the function Show[ ]:

Show[{ContourPlot3D[{(x^2) y==z,(x^2)+(y^2)==1},{x,-2,2},{y,-2,2},{z,-2,2},Mesh->False,PlotPoints->20,ContourStyle->{{RGBColor[0.501961, 1, 1],Specularity[1,100]},{Yellow,Opacity[0.4]}},Lighting->"Accent",BoundaryStyle->{1->None,2->None,{1,2}->{Red,Thick}}],
Graphics3D[{Arrowheads[.05],Arrow[Tube[{{0,0,0},{0,0,1.5}}]],
Arrow[Tube[{{0,0,0},{0,1.5,0}}]],
Arrow[Tube[{{0,0,0},{1.5,0,0}}]]}]}]

Finally, we do not want Axes, do not want it Boxed, and we want the Background to be Black. We set the opacity α of the cylinder to zero. With SphericalRegion -> True we don’t get jumping boundaries when we rotate the graphics with the mouse.

So the final code to make this graphic is just a few lines:

Show[{ContourPlot3D[{(x^2) y==c,(x^2)+(y^2)==1},{x,-2,2},{y,-2,2},{c,-2,2.},Mesh->False,Axes->False,Boxed->False,Background->Black,PlotPoints->20,SphericalRegion->True,ContourStyle->{{RGBColor[0.501961, 1, 1],Specularity[1,100]},{Yellow,Opacity[0]}},Lighting->"Accent",BoundaryStyle->{1->None,2->None,{1,2}->{Red,Thick}}],
Graphics3D[{Arrowheads[.05],Arrow[Tube[{{0,0,0},{0,0,1.5}}]],
Arrow[Tube[{{0,0,0},{0,1.5,0}}]],
Arrow[Tube[{{0,0,0},{1.5,0,0}}]]}]}]

Extensions

The possibilities in Mathematica are endless.

In the next example the surface is textured with the Texture of a cloud image 😀!
And the boundary contour is Pink and Dashed.

Show[{ContourPlot3D[{(x^2) y==c,(x^2)+(y^2)==1},{x,-2,2},{y,-2,2},{c,-2,2.},Mesh->False,Axes->False,Boxed->False,Background->Black,PlotPoints->20,SphericalRegion->True,Lighting->"Accent",ContourStyle->{Texture[],{Yellow,Opacity[0]}},BoundaryStyle->{1->None,2->None,{1,2}->Directive[Pink,Thick,Dashed]}],Graphics3D[{Arrowheads[.05],Arrow[Tube[{{0,0,0},{0,0,1.5}}]],
Arrow[Tube[{{0,0,0},{0,1.5,0}}]],
Arrow[Tube[{{0,0,0},{1.5,0,0}}]]}]}]

We can even go a step further, and modify the boundary Line segments into Tube segments.

We simply use the function ReplaceAll[ ] (with short notation /. ) to replace any occurrence of a Line[x_] function (where the underscore _ is ‘any element’, and x_ is ‘any element named x’) into a Tube[x, radius] function, where x is used for the element positions.

ContourPlot3D[{(x^2) y==c,(x^2)+(y^2)==1},{x,-2,2},{y,-2,2},{c,-2,2.},Mesh->False,Axes->False,Boxed->False,Background->Black,PlotPoints->20,SphericalRegion->True,ContourStyle->{{RGBColor[0.501961, 1, 1],Specularity[1,100]},{Yellow,Opacity[0.4]}},Lighting->"Accent",BoundaryStyle->{1->None,2->None,{1,2}->Red}]/. Line[x_]:>Tube[x,0.03]

We now Manipulate the tube radius and the opacity of the cylinder:

Manipulate[Show[{ContourPlot3D[{(x^2) y==c,(x^2)+(y^2)==1},{x,-2,2},{y,-2,2},{c,-2,2.},Mesh->False,Axes->False,Boxed->False,Background->Black,PlotPoints->20,ContourStyle->{{RGBColor[0.501961, 1, 1],Specularity[1,100]},{Yellow,Opacity[opacity]}},Lighting->"Accent",BoundaryStyle->{1->None,2->None,{1,2}->Red}]/. Line[x_]:>Tube[x,tuberadius],Graphics3D[{Arrowheads[.05],Arrow[Tube[{{0,0,0},{0,0,1.5}}]],
Arrow[Tube[{{0,0,0},{0,1.5,0}}]],
Arrow[Tube[{{0,0,0},{1.5,0,0}}]]}]}],{{tuberadius,0.1},0.01,0.2},{{opacity,0.5},0,1}]

It works beautifully! And the graph can smoothly be rotated, zoomed (hold the Ctrl key) and panned (hold the ⇧ key) with the mouse. Right-click gives save, print and reset options.

Intermezzo: Explanation of ReplaceAll

A bit explanation about the Replace operation:

/. means: ReplaceAll. An example:

{a,b,c,d,b}/.b->4

The underscore _ means ‘anything’. And x_ means: ‘anything, and we give it the name x’. It is used in functions:

f[x_]:=(x^2)

f[3]

f[monkey]

The := means ‘delayed assignment’, i.e. execute the function not at the definition time, but only then when it is called.

f1[x_]=Date[]

f2[x_]:=Date[]

Study the difference:

f1[x]

f2[x]

In the same way :> is a ‘RuleDelayed’. Alias: Esc:>Esc.
(PS: Esc is entered as \[EscapeKey ], without the space. See also: https://reference.wolfram.com/language/tutorial/KeyboardShortcutListing.html).

Generating a movie of the slowly rotating graph

For the movie, the frames are generated as a Table. We change the ViewPoint along a circle with radius 3, at level z=1. To speed up the computations, Parallelize takes care of automatically distributing the work over (in this case) 4 cores of the CPU of the Raspberry Pi 5. On your laptop you may be able to use more cores, and the CPU may be faster.

The SphericalRegion is defined as a sphere, with the origin of the CountourPlot in the center and with radius=2.2.

Timing[Parallelize[
movie=Table[Show[{ContourPlot3D[{x^2 y==c,x^2+y^2==1},{x,-2,2},{y,-2,2},{c,-2,2.},Mesh->False,Axes->False,Boxed->False,Background->Black,PlotPoints->20,SphericalRegion->True,ViewPoint->{3Cos[ϕ],3Sin[ϕ],1},Lighting->"Accent",ContourStyle->{Texture[],{Yellow,Opacity[0.5]}},BoundaryStyle->{1->None,2->None,{1,2}->Directive[Pink,Thick]}],Graphics3D[{Arrowheads[.05],Arrow[Tube[{{0,0,0},{0,0,1.5}}]],Arrow[Tube[{{0,0,0},{0,1.5,0}}]],Arrow[Tube[{{0,0,0},{1.5,0,0}}]]}]}],{ϕ,-π/2-0.1π,-π/2+0.1π,0.005π}];]]

ListAnimate[movie,AnimationDirection->ForwardBackward]

Of course, a finer time step gives a smoother movie.

Graphical User Interfaces (GUI)

A fancy GUI for a weather station

In Mathematica interesting and complex GUIs can be constructed. 
Here is an example of a weather station GUI, and of course there are a zillion of other possibilities.

First we collect weather data with Mathematica’s powerful functions to read external databases and websites:

getWeatherData[]:=Module[{temperature,pressure,pressureYesterday,pressureLinesMap,windRadar,humidity,sunrise,sunset,moonPhase,pressureHistory},
temperature=WeatherData["Eindhoven","Temperature"];
pressure=WeatherData["Eindhoven","Pressure"];
pressureYesterday=Normal[WeatherData["Eindhoven","Pressure",DatePlus[Today,-1]]][[1,2]];
pressureLinesMap=Import["https://www.knmi.nl/nederland-nu/weer/waarschuwingen-en-verwachtingen/weerkaarten","Images"][[4]];
windRadar=Import["https://www.weer.ws/windradar/","Images"][[4]];
humidity=WeatherData["Eindhoven","Humidity"];
sunrise=DateString[TimeZoneConvert[TimeObject[Sunrise[]]],{"Hour",":","Minute"}];
sunset=DateString[TimeZoneConvert[TimeObject[Sunset[]]],{"Hour",":","Minute"}];
moonPhase=Import["https://timesprayer.com/en/moon/netherlands-nl/amsterdam/","Images"][[3]];
pressureHistory=DateListPlot[[WeatherData["Eindhoven","Pressure",{DatePlus[Today,-14],Today}]], LabelingFunction->(Placed[Last@#1,Tooltip] &),
 Filling->Bottom,AspectRatio->0.3,ImageSize->300,PlotLabel->"Pressure last 14 days (mBar)"];
{temperature,pressure,pressureYesterday,pressureLinesMap,windRadar,humidity,sunrise,sunset,moonPhase,pressureHistory}];

Fetching all current weather data takes about a minute. This shows what is collected (give the following command a semicolon ; at the end, and no output will be shown).

See Mathematica busy scraping for you:

Timing[weatherData=getWeatherData[]]

The following Graphics canvas is Framed with a thick purple border with rounded corners, and filled with a collection of Inset functions, where we can precisely specify its content and location. The whole GUI is updated continuously, while the weather data are scraped every hour. The wind radar is shown as a movie.

task=RunScheduledTask[weatherData=getWeatherData[],3600]

Framed[Graphics[{Lighter[Yellow,0.85],Rectangle[{0,0},{192,108}],Inset[Style[DateString[{"DayName",", ","Day"," ","MonthName"," ","Year"}],30,Bold,Black,FontFamily->"EB Garamond 08"],{96,108}],Inset[Dynamic[Refresh[ClockGauge[Now,ImageSize->{100,100}],UpdateInterval->1]],{180,91}],Inset[Dynamic[Refresh[Style[DateString[TimeObject[]],45,Bold,Red,FontFamily->"Bitstream Vera Sans Mono"],UpdateInterval->1]],{114,90}],
Inset[Dynamic[Show[ThermometerGauge[weatherData[[1]],{-5,40},GaugeLabels -> "  Outside"],ImageSize->40]],{84,52}],
Inset[Show[ThermometerGauge[Quantity[20,"DegreesCelsius"],{0,40},GaugeLabels->"Inside"],ImageSize->40],{72,52}],
Inset[Dynamic[Show[weatherData[[4]],ImageSize->160]],{126,54}],
Inset[n=1;Dynamic[Pause[0.5];n++;Show[weatherData[[5]][[Mod[n,12]+1
]],ImageSize->200]],{23,57}],
Inset[Dynamic[AngularGauge[{weatherData[[2]],weatherData[[3]]},{950,1050},ScaleOrigin->{π,0},ScalePadding->.2,GaugeLabels->Row[{ Style["Yesterday",Bold,Blue],"  mBar  ",Style["Today",Bold,Red]}],PlotTheme->"Scientific",ImageSize->100]],{180,52}],
Inset[Dynamic[Show[AngularGauge[{40,100 weatherData[[6]]},{0,100},ScaleOrigin->{π,0},ScalePadding->0.15,PlotTheme->"Default",GaugeLabels->Row[{ Style["Inside",Bold,Red],"   Humidity   ",Style["Outside",Bold,Blue]}]],ImageSize->80]],{60,1}],
Inset[Dynamic[Style["Sunrise: "<>weatherData[[7]]<>"\nSunset : "<>weatherData[[8]],Black,Bold,10]],{0,0}],
Inset[Dynamic[Show[weatherData[[9]],ImageSize->60]],{30,0}],
Inset[Dynamic[weatherData[[10]]],{140,10}]
},Frame->False,ImageSize->3{192,108}],RoundingRadius->50,Background->Lighter[Yellow,0.85],FrameStyle->Directive[Purple,Thickness[15]]]

Be sure to remove the scheduled task when it’s no longer needed to free up system resources.

RemoveScheduledTask[task]

Physics

Mesmerizing bouncing balls - the visible physics of collisions

The source idea comes from a wonderful series of explaining pages by Matthias Müller:

Ten Minute Physics: https://matthias-research.github.io/pages/tenMinutePhysics/index.html

See: 03 - Ball collision handling in 2d, YouTube: https://www.youtube.com/watch?v=ThhdlMbGT5g

I asked the Notebook Assistant to convert the HTML code into Wolfram Language code, and the result was pleasing! The names of the variables are such that the code is more or less self-explaining.

Initialization

ClearAll[simulation,drawScene,handleCollisions];

Parameters

timeStep=1/60;
gravitationalAcceleration={0,2}; (*{x,y}*)
worldSize={3,3}; (*Simulation area width and height*)
restitution=0.95;
numberOfBalls=20;

Ball structure: {position,velocity,radius,mass}

For every ball, the initial position, velocity, radius and mass is set randomly, each within a specific range:

ballsInit:=Table[{
RandomReal/@worldSize,(*Initial position*)
RandomReal[{-1,1},2],(*Initial velocity*)
RandomReal[{0.05,0.1}],(*Radius*)
RandomReal[{0.1,1}] (*Mass*)
},{numberOfBalls}];

balls=ballsInit;

Function to simulate the system

In the simulation step the positions and velocities of the balls are updated. Note that when a ball comes nearer to another ball than the sum of their radius values, the velocity vector gets a negative sign: the ball goes further in the other way after a bump.

simulation[]:=(
balls=Map[Function[ball,Module[{pos,vel,radius,mass,newVel},{pos,vel,radius,mass}=ball;
(*Update velocity with gravity*)newVel=vel-gravitationalAcceleration timeStep;
(*Update position*)pos=pos+newVel timeStep;
(*Handle collisions with walls*)If[pos[[1]]-radius<0||pos[[1]]+radius>worldSize[[1]],newVel[[1]]=-newVel[[1]];
pos[[1]]=Clip[pos[[1]],{radius,worldSize[[1]]-radius}]];
If[pos[[2]]-radius<0||pos[[2]]+radius>worldSize[[2]],newVel[[2]]=-newVel[[2]];
pos[[2]]=Clip[pos[[2]],{radius,worldSize[[2]]-radius}]];
{pos,newVel,radius,mass}]],balls];
(*Handle ball-ball collisions*)
handleCollisions[];);

Function to handle ball-ball collisions

handleCollisions[]:=(Do[Module[{pos1,vel1,r1,m1,pos2,vel2,r2,m2,dir,dist,overlap,v1,v2},{pos1,vel1,r1,m1}=balls[[i]];
{pos2,vel2,r2,m2}=balls[[j]];
dir=pos2-pos1;
dist=Norm[dir];
If[dist>0&&dist<r1+r2,dir=Normalize[dir];
overlap=(r1+r2-dist)/2;
pos1-=overlap dir;
pos2+=overlap dir;
v1=Dot[vel1,dir];
v2=Dot[vel2,dir];
{v1,v2}={(m1 v1+m2 v2-m2 (v1-v2) restitution)/(m1+m2),(m1 v1+m2 v2-m1 (v2-v1) restitution)/(m1+m2)};
vel1+=(v1-Dot[vel1,dir]) dir;
vel2+=(v2-Dot[vel2,dir]) dir;
balls[[i]]={pos1,vel1,r1,m1};
balls[[j]]={pos2,vel2,r2,m2};]],{i,Length[balls]},{j,i+1,Length[balls]}];);

Function to draw the scene

drawScene[]:=Graphics[{EdgeForm[None],Table[{Purple,Disk[balls[[i,1]],balls[[i,3]]]},{i,Length[balls]}]},PlotRange->{{0,worldSize[[1]]},{0,worldSize[[2]]}},AspectRatio->Automatic];

Interactive Simulation

Here is the simulation. It runs remarkably well on our tiny computer. The ‘restitution‘ slider determines the amount of energy that is transferred between a bounce.

DynamicModule[{running=True},Panel[Column[{Row[{Button["Start",balls=ballsInit;running=True;],Button["Stop",running=False;]," Restitution: ",Slider[Dynamic[restitution],{0,1},ImageSize->50]}],Dynamic[If[running,simulation[];
Pause[0]];
drawScene[]]}]]]

Impressive commands on a tiny computer

Fluid flow simulation with Navier-Stokes equations modeling

Source: https://reference.wolfram.com/language/ref/StreamPlot.html

The Navier-Stokes equations are a set of partial differential equations that describe the motion of fluid substances such as liquids and gases. These equations account for various physical phenomena, including the effects of viscosity, pressure, and external forces on the fluid flow.

The general form of the Navier-Stokes equations can be expressed for an incompressible fluid as:

1. Momentum Conservation (Vector form):

Where:

(OverVector[u])={u,v} is the velocity vector field of the fluid.

is the pressure field.

is the density of the fluid.

is the dynamic viscosity.

represents external forces (e.g., gravity).

2. Mass Conservation (Continuity equation for incompressible flow):

(OverVector[∇])·(OverVector[u])=0

These equations are fundamental in fluid mechanics and are used to model weather patterns, ocean currents, water flow in pipes, air flow around wings, and many other applications.

For an explanation of the equations see (for example) this video: https://www.youtube.com/watch?v=Ra7aQlenTb8.

Solving the Navier-Stokes equations can be challenging due to their non-linearity and the complexity of the boundary conditions in practical scenarios.

It is easy in Mathematica:

Specify the geometry for a fluid flow:

Ω=RegionUnion[Rectangle[{0,0.01},{0.02,0.02}],Rectangle[{0.02,0},{0.08,0.03}],Rectangle[{0.08,0},{0.1,0.01}],Rectangle[{0.08,.02},{0.1,0.03}],Rectangle[{0.04,0.03},{0.06,0.04}]];
RegionPlot[Ω,AspectRatio->Automatic]

Specify the boundary conditions (BC) for a laminar flow that enters at the left-hand channel and exits through the two right-hand channels:

inflowBC=DirichletCondition[{u[x,y]==4000(y-0.01)(0.02-y),v[x,y]==0},x==0];outflowBC=DirichletCondition[p[x,y]==0,x==0.1];wallBC=DirichletCondition[{u[x,y]==0,v[x,y]==0},0<x<0.1];bcs={inflowBC,outflowBC,wallBC};

Define the Navier–Stokes equations for u, v and mass conservation (note the beautiful mathematical notation in the input code):

navierstokes={
(Div[(-μ (Grad[u[x,y], {x,y}])), {x,y}])+ρ {u[x,y],v[x,y]}.(Grad[u[x,y], {x,y}])+Derivative[1,0][p][x,y]==0,(Div[(-μ (Grad[v[x,y], {x,y}])), {x,y}])+ρ {u[x,y],v[x,y]}.(Grad[v[x,y], {x,y}])+Derivative[0,1][p][x,y]==0,
Derivative[1,0][u][x,y]+Derivative[0,1][v][x,y]==0
}/. {ρ->1,μ->(2/(10^5))};

Use the finite element method to solve for the steady flow velocities and pressure.
It is solved in less than a minute on a Raspberry Pi 5:

Timing[{uVel,vVel,pressure}=NDSolveValue[{navierstokes,bcs},{u,v,p},{x,y}∈Ω,Method->{"FiniteElement","InterpolationOrder"->{u->2,v->2,p->1},"MeshOptions"->{"MaxCellMeasure"->0.00000025}}];]

Use (RefLink[DensityPlot, paclet:ref/DensityPlot]) to plot the pressure and (RefLink[StreamPlot, paclet:ref/StreamPlot]) to plot the flow:

Show[
DensityPlot[pressure[x,y],{x,y}∈uVel["ElementMesh"],PlotRange->All,Frame->None,Axes->None,AspectRatio->Automatic,ColorFunction->"DeepSeaColors",PlotLegends->Placed[BarLegend[Automatic,LegendLabel->Placed["pressure",Above]],Above]],StreamPlot[{uVel[x,y],vVel[x,y]},{x,y}∈Ω,AspectRatio->Automatic,StreamColorFunction->"SolarColors",PlotLegends->BarLegend[Automatic,LegendLabel->Placed["speed",Below]],StreamPoints->Fine],ImageSize->Full]

Pisot numbers: amazing and remarkable

Mathematica is very precise:

N[π,50000]

N[(((((((2/(27+3 (Sqrt[69])))))^(1/3))+(1/3) ((((1/2) (27+3 (Sqrt[69]))))^(1/3))))^30000),10000]

Note the many consecutive zero's in this remarkable Pisot number.

How to Parallelize my code with a single command?

To see the 4 CPUs of our Raspberry Pi 5 working simultaneously when we run Mathematica’s fully automatic Parallelize[ ] function, we install the Gnome system monitor. Type > as the first character to start a cell that evaluates in another language (click on the icon to the left of the cell). Here we run a Linux shell command to install the program:

sudo apt install gnome-system-monitor -y

Start the Gnome system monitor on your RPI5, and select in its top-bar ‘Resources’.:

gnome-system-monitor

As this cell keeps running, abort with Alt. (alt-key + point), or use the Abort key in the upper left of this notebook window (with the centered black dot).

In Windows: Open Task Manager (CtrlAlt⌫) ▶ Performance ▶ Open Resource Monitor to visually check the CPU activity.

Now we start a parallel calculation: here a test if a high exponent of two minus 1 is a prime. The Q in PrimeQ means ‘question mark’. Note the start of the 4 CPUs when you start the execution:

res=Parallelize[Map[PrimeQ[(2^(#))-1]&,Range[9601,12000]]];

Short[res,10]

A long list (2400 items) of False outcomes, and maybe there are a few True findings. 
Find their position with the following command:
(as we started at 9601, and the function Position[ ] starts at 1, we need to add 9600 to the positions):

Position[res,True]+9600

Just check:

PrimeQ[(2^11213)-1]

It must be a big number:

(2^11213)

Watch the results appear as they are found, now we take the range 1000-5000:

SetSharedVariable[primes]

primes={};
Monitor[Parallelize[Scan[If[PrimeQ[(2^(#))-1],AppendTo[primes,#]]&,Range[1000,5000]],Method->"FinestGrained"],primes];
primes

If you want to find out more about the organization of the parallel kernels on your system:
Click in the front end menu bar Evaluation ▶ Parallel Kernel Configuration

Conclusion: you can now do EVERYTHING in 21st century coding

STEM students: get your RPI5 or a student license (cheaper than Netflix)

Parents, retired, schools: give / use a RPI5 for fun and learning

Current LLMs today typically learn Python, C++ and Java. It would be timely to add millions of Mathematica users! Spread the word, start your YouTube channel, set up your school team!

Check the resources:

Documentation (best in the world): fun to keep learning! https://reference.wolfram.com/language/

Public resources: www.wolfram.com/resources/

Thousands of cool demonstrations: https://demonstrations.wolfram.com/

The notebook of this lecture gives much more details

A gentle introduction to Mathematica (in 1 hour) and many applications: www.21stcenturycoding.net,
see also www.youtube.com/@21stcenturycoding

A personal note: I have been a STEM student my whole life

I feel privileged I could work with Mathematica notebooks for 35 years, 
brain-inspired computing for medical imaging (big data), it made my career

Spread the word of these STEM delights
in the Raspberry Pi communities,
social media, schools and beyond!

Appendix: Radial-Basis-Transform

Radial thin-plate-spline functions:

$HistoryLength=0;
Clear[ThinPlateSpline];
ThinPlateSpline[1] = Power[#, 3/2] &;
ThinPlateSpline[2] = (# Log[# + $MinMachineNumber] / 2)&;
ThinPlateSpline[3] = Sqrt;
ThinPlateSpline[_Integer?Positive] = Identity;

A fast internal utility function to calculate the distance matrix between two sets of vectors:

ClearAll[InternalDistanceMatrix];

InternalDistanceMatrix := 
	InternalDistanceMatrix = 
		Replace[
			LibraryFunctionLoad[
				FindLibrary["NumericArrayUtilities"],
				"distance_matrix_MTensor2Arg", 
				{{Real, 2, "Constant"}, {Real, 2, "Constant"}, Integer, True | False},
				{Real, 2}
			],
			Except[_LibraryFunction] :> 
			GeneralUtilities`Panic["RadialBasisTransformError", "Couldn't load \"distance_matrix_MTensor2Arg\" from NumericArrayUtilities."]
		];

Implementation of the RadialBasisTransform:

Parameters A and (Subscript[w, i]) are estimated by solving the linear system (({{B-λ I, Q}, {(Q^T), 0}}))(({{w}, {A}}))=(({{y}, {0}})), where B=((Subscript[b, ij]))=ϕ((Subscript[x, i])-(Subscript[x, j])), λ is the smoothness parameter, and Q=(({{1, (Subscript[x, 1,1]), (Subscript[x, 1,2]), …, (Subscript[x, 1,d])}, {1, (Subscript[x, 2,1]), (Subscript[x, 2,2]), …, (Subscript[x, 2,d])}, {⋮, ⋮, ⋮, ⋱, ⋮}, {1, (Subscript[x, n,1]), (Subscript[x, n,2]), …, (Subscript[x, n,d])}})) with d being the number of dimensions and n the number of basis points.

ClearAll[RadialBasisTransform];

RadialBasisTransform::usage = 
"RadialBasisTransform[{(Subscript[p, 1])->(Subscript[q, 1]),…},rbf] 
gives a continuous transformation mapping point (Subscript[p, i]) to (Subscript[q, i]) using the radial basis function rbf.\n RadialBasisTransform[{(Subscript[p, 1]),…}->{(Subscript[q, 1]),…},rbf]
gives a transformation mapping sources {(Subscript[p, 1]),…} to targets {(Subscript[q, 1]),…}.\n RadialBasisTransform[{(Subscript[p, 1]),…}->{(Subscript[q, 1]),…},rbf,λ]
returns for λ>0 a regularized approximation rather than an interpolation.";

Options[RadialBasisTransform] = {};

RadialBasisTransform[
	map:Except[_String], 
	radialbasis:(_String|{_String, _Integer}):"ThinPlateSpline", 
	λ_:Automatic,
	opts:OptionsPattern[RadialBasisTransform]
] := 
Catch[
	Module[
		{sources, targets, eps, df = "SquaredEuclideanDistance", lambda, n, sdim, tdim, rbf, Q, B, M, Wt, a0, At},
		
		(* anchor point map *)
		Switch[
			map,
			_List -> _List, {sources, targets} = List @@ map,
			{__Rule}, {sources, targets} = Transpose[List @@@ map],
			_, Message[RadialBasisTransform::anchortarget]; Throw[$Failed, RadialBasisTransform]
		];
		If[
			Length[sources]=!=Length[targets],
			Message[RadialBasisTransform::anchordims]; Throw[$Failed, RadialBasisTransform]
		];
		If[
			Not[ArrayQ[sources, 1|2, Internal`RealValuedNumericQ] && ArrayQ[targets, 1|2, Internal`RealValuedNumericQ]],
			Message[RadialBasisTransform::anchortarget]; Throw[$Failed, RadialBasisTransform]
		];
		{n, sdim} = PadRight[Dimensions[sources], 2, 0];
		tdim = Last@PadRight[Dimensions[targets], 2, 0];
		If[sdim===0, sources=Thread[{sources}]];
		If[tdim===0, targets=Thread[{targets}]];
		eps = 1/Norm[Flatten[Map[Differences, Map[MinMax, Transpose[sources]]]]/(n^(1/Max[sdim, 1])+2)];
		
		(* radial basis function *)
		Switch[
			radialbasis,
			"ThinPlateSpline"|Automatic,
				rbf = ThinPlateSpline[Max[sdim, 1]],
			{"ThinPlateSpline"|Automatic, _?Positive},
				rbf = ThinPlateSpline[Max[sdim, 1]]; eps = N@Last[radialbasis],
			{"PolyharmonicSpline", _Integer?(Positive[#] && OddQ[#]&)},
				rbf = #^(radialbasis[[2]]/2)&,
			{"PolyharmonicSpline", _Integer?(Positive[#] && EvenQ[#]&)},
				rbf = #^(radialbasis[[2]]/2) Log[#]/2&,
			{"PolyharmonicSpline", _Integer?(Positive[#] && OddQ[#]&), _?Positive},
				rbf = #^(radialbasis[[2]]/2)&; eps = N@Last[radialbasis],
			{"PolyharmonicSpline", _Integer?(Positive[#] && EvenQ[#]&), _?Positive},
				rbf = #^(radialbasis[[2]]/2) Log[#]/2&; eps = N@Last[radialbasis],
			"Gaussian",
				rbf = Exp[-#]&,
			{"Gaussian", _?Positive},
				rbf = Exp[-#]&; eps = N@Last[radialbasis],
			"MultiQuadratic",
				rbf = Sqrt[1 + #]&,
			{"MultiQuadratic", _?Positive},
				rbf = Sqrt[1 + #]&; eps = N@Last[radialbasis],
			"InverseQuadratic",
				rbf = 1/(1 + #)&,
			{"InverseQuadratic", _?Positive},
				rbf = 1/(1 + #)&; eps = N@Last[radialbasis],
			"InverseMultiQuadratic",
				rbf = 1/Sqrt[1 + #]&,
			{"InverseMultiQuadratic", _?Positive},
				rbf = 1/Sqrt[1 + #]&; eps = N@Last[radialbasis],
			{"Bessel", _Integer?(Positive[#]&)},
				rbf = (BesselJ[radialbasis[[2]]/2 - 1, #] #^-(radialbasis[[2]]/2 - 1))&,
			{"Bessel", _Integer?(Positive[#]&), _?Positive},
				rbf = (BesselJ[radialbasis[[2]]/2 - 1, #] #^-(radialbasis[[2]]/2 - 1))&; eps = N@Last[radialbasis],
			_,
				Message[RadialBasisTransform::nrbf, radialbasis]; Throw[$Failed, RadialBasisTransform]
		];
		
		(* lambda *)
		Switch[
			λ,
			_?NonNegative, lambda = N[λ],
			Automatic, lambda = Automatic,
			_, Message[RadialBasisTransform::nlambda, λ]; Throw[$Failed, RadialBasisTransform]
		];
		
		If[
			Max[1, tdim] + 1 === n,
			(* create linear transformation *)
			Q = PadLeft[N@sources, {n, Max[sdim, 1] + 1}, 1.`];
			{{a0}, At} = Internal`PartitionRagged[PseudoInverse[Q].targets, {1, sdim}];
			RadialBasisTransform[
				"RadialBasisFunction", 
				radialbasis, 
				sdim->tdim,
				sources,
				Which[
					tdim === 0,
					With[
						{a$0=First[a0], A$t=Flatten[At]},
						a$0+#.A$t&
					],
					tdim === 1,
					With[
						{a$0=First[a0], A$t=Flatten[At]},
						Transpose[{a$0+#.A$t}]&
					],
					tdim >= 2,
					With[
						{a$0=a0, A$t=At},
						Transpose[a$0+Transpose[#.A$t]]&
					]
				]
			],
			(* create nonlinear transformation *)
			Q = PadLeft[N@sources, {n, Max[sdim, 1]+1}, 1.`];
			B = rbf[InternalDistanceMatrix[#, #, Replace[df, {"EuclideanDistance"->3, "SquaredEuclideanDistance"->4}], True]&[eps * N@sources]];
			(* determine Automatic value of lambda here *)
			If[lambda === Automatic, lambda = 0.`];
			M = ArrayFlatten[{{B - lambda IdentityMatrix[n], Q}, {Transpose[Q], ConstantArray[0.`, {Max[sdim, 1], Max[sdim, 1]}+1]}}];
			Off[LinearSolve::krymit];
			{Wt, {a0}, At} = Internal`PartitionRagged[LinearSolve[M, PadRight[N@targets, {n+Max[sdim, 1]+1, Max[tdim, 1]}, 0.`], Method-> "Krylov"], {n, 1, Max[sdim, 1]}];
			On[LinearSolve::krymit];
			RadialBasisTransform[
				"RadialBasisFunction", 
				radialbasis, 
				sdim->tdim,
				sources,
				Which[
					tdim === 0,
					With[
						{a$0=First[a0], A$t=Flatten[At], rbf$=rbf, sources$=eps*sources, W$t=Flatten[Wt], eps$=eps, df$=Replace[df, {"EuclideanDistance"->3, "SquaredEuclideanDistance"->4}]},
						a$0+#.A$t+rbf$[InternalDistanceMatrix[eps$*#, sources$, df$, True]].W$t&
					],
					tdim === 1,
					With[
						{a$0=First[a0], A$t=Flatten[At], rbf$=rbf, sources$=eps*sources, W$t=Flatten[Wt], eps$=eps, df$=Replace[df, {"EuclideanDistance"->3, "SquaredEuclideanDistance"->4}]},
						Transpose[{a$0+#.A$t+rbf$[InternalDistanceMatrix[eps$*#, sources$, df$, True]].W$t}]&
					],
					tdim >= 2,
					With[
						{a$0=a0, A$t=At, rbf$=rbf, sources$=eps*sources, W$t=Wt, eps$=eps, df$=Replace[df, {"EuclideanDistance"->3, "SquaredEuclideanDistance"->4}]},
						Transpose[a$0+Transpose[#.A$t+rbf$[InternalDistanceMatrix[eps$*#, sources$, df$, True]].W$t]]&
					]
				]
			]
		]
	],
	RadialBasisTransform
]

Patterns for RadialBasisTransform execution:

RadialBasisTransform /: RadialBasisTransform["RadialBasisFunction", rbf_, dims_Rule, sources_List, f_Function][prop_String]:=
Switch[
	prop,
	"InterpolationMethod",
		"RadialBasisFunction",
	"InputDimension",
		First[dims],
	"OutputDimension",
		Last[dims],
	"Coordinates",
		If[First[dims]===0,Flatten[#],#]&[sources],
	"Values"|"ValuesOnCoordinates",
		RadialBasisTransform["RadialBasisFunction", dims,f][
			If[First[dims]===0,Flatten[#],#]&[sources]
		],
	"RadialBasis", rbf,
	_,
		Undefined
]

RadialBasisTransform /: RadialBasisTransform["RadialBasisFunction", rbf_, 0 -> 0, sources_List, f_Function][arg_?NumericQ] := First[f[{{arg}}]]
RadialBasisTransform /: RadialBasisTransform["RadialBasisFunction", rbf_, 0 -> 0, sources_List, f_Function][args:{__?NumericQ}] := f[Transpose@{args}]
RadialBasisTransform /: RadialBasisTransform["RadialBasisFunction", rbf_, 0 -> _Integer?Positive, sources_List, f_Function][arg_?NumericQ] := First[f[{{arg}}]]
RadialBasisTransform /: RadialBasisTransform["RadialBasisFunction", rbf_, 0 -> _Integer?Positive, sources_List, f_Function][args:{__?NumericQ}] := f[Transpose@{args}]
RadialBasisTransform /: RadialBasisTransform["RadialBasisFunction", rbf_, _Integer?Positive -> 0, sources_List, f_Function][arg:{__?NumericQ}] := First[f[{arg}]]
RadialBasisTransform /: RadialBasisTransform["RadialBasisFunction", rbf_, _Integer?Positive -> 0, sources_List, f_Function][args:{{__?NumericQ}..}] := f[args]
RadialBasisTransform /: RadialBasisTransform["RadialBasisFunction", rbf_, _Integer?Positive -> _Integer?Positive, sources_List, f_Function][arg:{__?NumericQ}] := First[f[{arg}]]
RadialBasisTransform /: RadialBasisTransform["RadialBasisFunction", rbf_, _Integer?Positive -> _Integer?Positive, sources_List, f_Function][args:{{__?NumericQ}..}] := f[args]

Typesetting of RadialBasisTransform:

RadialBasisTransform /: MakeBoxes[a_RadialBasisTransform, fmt_] := BoxForm`MakeRowBoxes[a, fmt] /; Refresh[BoxForm`UseTextFormattingQ, None];

RadialBasisTransform /: MakeBoxes[itrafo_RadialBasisTransform, fmt_] /; And[BoxForm`UseIcons] := 
Module[
	{icon, idim, odim, anchors, rbf,alwaysGrid},
	idim = Replace[itrafo @ "InputDimension", 0->"scalar"];
	odim = Replace[itrafo @ "OutputDimension", 0->"scalar"];
	rbf= itrafo @ "RadialBasis";
	icon = BoxForm`GenericIcon[InterpolatingFunction];
	anchors = Length[itrafo @ "Coordinates"];
	alwaysGrid = {
		BoxForm`SummaryItem @ {"Radial basis: ", rbf},
		BoxForm`SummaryItem @ {"Dimensions: ", idim -> odim},
		BoxForm`SummaryItem @ {"# of anchors: ", anchors}
	};
	BoxForm`ArrangeSummaryBox[RadialBasisTransform, itrafo, icon, alwaysGrid, {}, fmt]
];

Boolean token to speed up ImageTransformation:

RadialBasisTransform /: Image`SpatialOperationsDump`ListableQ[_RadialBasisTransform] = True;

NotebookFind[InputNotebook[],"Output",All,CellStyle];
FrontEndTokenExecute[InputNotebook[],"SelectionConvert","BitmapConditional"];